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We study the immersion of a ferromagnetic nanowire within a nematic liquid crystal using a lattice 
Boltzmann algorithm to solve the full three-dimensional equations of hydrodynamics. We present 
an algorithm for including a moving boundary, to simulate a nanowire, in a lattice Boltzmann 
simulation. The nematic imposes a torque on a wire that increases linearly with the angle between 
the wire and the equilibrium direction of the director field. By rotation of these nanowires, one can 
determine the elastic constants of the nematic. 
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I. INTRODUCTION 

In a pioneering theoretical study in 1973 de Gennes 
and Brochard looked at the effect of dispersing long rod- 
like particles in a liquid crystal pj. Despite interesting 
predictions on ferronematics and fcrrocholesterics, little 
experimental work was done at the time. Very recently, 
however, there has been a resurgence of interest in this 
subject motivated by the advent of nanofabrication tech- 
niques to make such particles 0, El 0, IE • Early results 
are already very encouraging. For example, solutions of 
just 0.1 — 0.2wt% of carbon nanotubes in E7 have shown a 
dramatic increase (~ 30°) in the nematic-isotropic tran- 
sition temperature 0, IU • 

Recent research on the inclusion of spherical particles 
in a nematic have shown that unique interactions be- 
tween particles which can lead to self-assembly [ol ITnllllj . 
There has also been considerable interest in generalizing 
micro-rheological tools, such as examining the fluctua- 
tions of microspheres to measure local viscosities |l2ill3j . 
However, there are complications in using spheres for 
micro-rheology in an anisotropic fluid like a liquid crystal 
that are not present in isotropic fluids. First, boundary 
conditions at the surface of the sphere typically induce 
defects either on or close to the surface of the sphere 
[Hi El Il5l El • These significantly modify the director 
field of the liquid crystal thus making the sphere far from 
a passive probe. Second, these spheres are typically ma- 
nipulated via laser traps However, the elec- 
tric fields required to trap a sphere also strongly affect 
the alignment of the liquid crystal |20| again making the 
system non-passive. 

Recent experiments on nickel nanowires in a liquid 
crystal have shown that they can be manipulated with 
magnetic fields and that the torque exerted on the wire 
by the liquid crystal can be directly measured 0, |2^| . 
The magnetic fields required to manipulate these ferro- 
magnetic wires are typically very small, ~ 5 Gauss [2l| . 
much less than that required to align the liquid crystal 
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FIG. 1: Schematic of a sample cell in our simulation. The 
dimensions are 1.25/im x 1.25/xm x 1.25/im. Two rigid sur- 
faces are located at z — and 2 = L and there is periodic 
boundaries in the xy— plane. In the middle of the cell is a 
ferromagnetic nanowire 625nm in length and approximately 
88nm in diameter. The top front quarter has been removed 
for clarity. 



(~ 3500 Gauss for the Frederik's transition in 5CB in a 
cell with thickness 10/xm|2J|). Further, when the wire 
is aligned with the director field, it causes no distortion. 
Thus, the use of such an anisotropic particle is ideal for 
doing micro-rheology in an anisotropic fluid. By manip- 
ulating the wires immersed in the nematic it should be 
possible to infer the local viscous and elastic properties 
of the liquid crystal. 

Simulations of colloids in liquid crystals have also 
recently been performed, but most have focused on 
homeotropic anchoring to the colloid. Monte Carlo sim- 
ulations were performed on rods of infinite length [24| . 
Also, there is work on soft two-dimensional cylinders 
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using lattice Boltzmann [25j and FEM simulations on 
hard cylinders |26f. All these works used homeotropic 
anchoring which leads to saturn ring-like defects along 
and around the cylinder. In this paper we focus on the 
cylinder, or wire, with boundary conditions where the 
nematic aligns with the long axis of the wire. The lack 
of defect formation from this configuration makes it ideal 
for microrheological applications as it can now measure 
the linear response of an unperturbed background. 

Analytic calculations of the response of a liquid crys- 
tal to an immersed particle can only be carried out by 
imposing rather unrealistic assumptions. This makes the 
use of simulations necessary to understand precisely what 
is measured in experiments. Lattice Boltzmann schemes 
have proved to be very successful in simulations of com- 
plex fluids 27], including liquid crystals |2l E^, We 
will examine a system similar to the experimental setup 
described in [2l|, as shown in Figure ^ We have a simu- 
lation cell with rigid surfaces for boundaries at z = and 
z = L, where L is the height of the cell. The surfaces im- 
pose a preferred aligned direction for the nematic liquid 
crystal. The xy— plane has periodic boundaries in both 
directions. At the center of the cell, a small ferromagnetic 
nickel wire will be immersed. With the system in equi- 
librium, we can rotate the magnetic field in the z plane 
to rotate the wire within the box. Using the response 
function of the imposed torque, we are able to measure 
the elastic properties of the nematic liquid crystal. 



II. EQUATIONS OF MOTION 

Liquid crystals can be characterized by the second 
moment of the orientation of the constituent molecules 
p, in terms of a local tensor order parameter Q where 
Qa/3 = (PaPp — \5 a p)- We use angular brackets to de- 
note a coarse-grained thermal average. Greek indices will 
be used to represent Cartesian components of vectors and 
tensors and the usual summation over repeated indices 
is assumed. Q is a traceless symmetric tensor. We label 
its largest eigenvalue |g so that < q < 1. q describes 
the magnitude of order along Q's principle eigenvector 
n, referred to as the director. 

The nematic order parameter evolves according to the 
convective-diffusive equation |30| 



(d t + u- V)Q - S(W, Q) = D r ll 



(1) 



where D r is a collective rotational diffusion constant. 
The order parameter distribution can be both rotated 
and stretched by flow gradients described by 

S(W,Q) = (xA + n)(Q + I/3) + (Q + I/3)( X A-n) 
-2 X (Q + I/3)Tr(A(Q + I/3)) (2) 

where A = (W + W T ) /2 and Q = (W - W T )/2 arc the 
symmetric part and the anti-symmetric part respectively 
of the velocity gradient tensor W a p = dpu a . x is a con- 
stant which depends on the molecular details of a given 
liquid crvstal|31|. 



The term on the right hand side of Eqn. Q describes 
the relaxation of the order parameter towards the min- 
imum of the free energy. The molecular field H that 
provides the driving motion is related to the variational 
derivative of the free energy by 



H=-— + ^Tr— 
3 5Q 



(3) 



The symmetry and zero trace of Q and H are exploited 
for simplification. 

The equilibrium properties of a liquid crystal can be 
described in terms of the Landau-de Gennes free energy 

m 

T= [ dV{F B + F E } + [ dSF s . (4) 

JV JdV 

Fb describes the bulk free energy 
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Fb — ttO- - T;)Qa0 ^—QupQpfQfa 



(•5) 



The bulk free energy describes a liquid crystal with a first 
order phase transition at 7 = 2.7 from the isotropic to 
nematic phase |22j. When 7 > 3 the isotropic phase is 
completely unstable. Unless noted, we use 7 = 3.2 and 
A = 1.5 x 10 5 pN/(fJ,m) 2 which leads to q = 0.555 in 
the unstressed equilibrium uniaxial case. We also used 
higher values of q and found no qualitative difference. 
The free energy also contains an elastic term, 



Li 



^-QaflidaQj^idpQye), (6) 

where the L's are the material specific elastic constants. 
It is straightforward to map L's to the Frank elastic con- 
stants (the splay elastic constant K\, the twist elastic 
constant K 2 , and the bend elastic constant if3)H3- The 
one elastic constant approximation K\ — K 2 = K3 = K 
can be easily achieved by setting Li > 0, L 2 = L3 = 
|33|. One advantage of simulations is the ease of switch- 
ing between such an approximation commonly used in 
analytical work and the full equations necessary to com- 
pare to experiment. 

At the surfaces of the system we assume a pinning 
potential 



Fs = \a s {Q aP -Q° a p) 



where Q° has the form 



(7) 



(8) 



and q° is set to the equilibrium bulk value q. This corre- 
sponds to a preferred direction n° for the director at the 
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surface. Typically we set a s = 6.33 x 10 4 pN/ /j,m, which 
corresponds to the strong pinning limit where Q ~ Q a at 
the boundaries. 

Inclusion of the wire in the system will be approxi- 
mated as a moving boundary. We will use a pining po- 
tential similar to (JJJ) for pining the director field of the 
liquid crystal to the long axis of the wire. This will be 
discussed in greater detail below and in Section lill Bl 

The fluid momentum obeys the continuity 



d t p + d a pu a = 
and the Navier-Stokes equation 



(9) 



p(d t + Up8p)u a = dpT a p + dpO a p 

+r)dp{(l - 3d p P )d 1 u~ f 5 al 3 + d a up + dpu a )(10) 

where p is the fluid density and rj — prf /3 is an isotropic 
viscosity. The form of this equation is not dissimilar to 
that for a simple fluid (anisotropic viscosities arise due 
to terms in a a p. However the details of the stress ten- 
sor reflect the additional complications of liquid crystal 
hydrodynamics. There is a symmetric contribution 



O a p 
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2£(<2a/3 + -^S a p)Q 7e H 7e - dpQ lv — q — (11) 



and an antisymmetric contribution 



(12) 



These additional terms can be mapped onto the Erickson- 
Leslie equations to give the Leslie coefficients 29J. The 
isostatic pressure P is constant in the simulations to 
a very good approximation, consistent with the liquid 
crystal being incompressible. 

At the surface of the wire we assume a molecular sur- 
face field of a form similar to Eqn J7J) where 



H w = w s (Q - Q w ), 



(13) 



and Q™a = q w (m a nip — S a p/3). The director m of Q u 
is fixed along the long axis of the wire and q w is set to 
the equilibrium bulk value 



(14) 



For the uniaxial case, where Q a p = q \n a np ^ 

the torque on the liquid crystal due to the presence of 
the wire is easily shown to be [2j| 

r = mxh"\ (15) 

where m is the unit direction vector of the wire, and 



q w {mpH^ +m a H^) 



In our simulations, the torque is determined in system 
variables by 

r Q = 2[Q a pH^ — Q a7 H™p + Hp^(Q(3p — Q 77 ) + 

Qp^H^-H^)], (17) 

which can be easily shown to be equivalent to Eqn. 115|) 
in the uniaxial case. 

The motion of the wire is performing a balancing act 
between forces. Its motion is described by 



A(m x rh) = —(T - (T ■ m)m) + /x x B 



(18) 



where A is a rotational viscosity and fi x B is the torque 
on the wire due it its magnetic moment \x in the applied 
magnetic field. The (r • m)m term is zero but it's pres- 
ence enhances the numerical stability of the scheme. In 
principle, A should depend on the orientation of the wire 
|22| . but in this study we are interested in the steady 
state solution. As such, this term only serves to move 
the wire between steady states here. 



III. NUMERICAL IMPLEMENTATION 

A. The Lattice Boltzmann Algorithm for Liquid 
Crystal Hydrodynamics 

We summarize a lattice Boltzmann scheme which 
reproduces Eqns. @, and l|lUfl to second order. 

The reader only interested in the results can safely 
skip over this section. We defined two distribution 
functions /i(x) and G^(x) where i labels the lattice 
directions from site x. In two dimensions, we choose 
a nine-velocity model on a square lattice with ve- 
locity vectors ~t t = (±1, 0), (0, ±1), (±1, ±1), (0, 0). 
For three-dimensional systems, a 15-velocity model 
on a cubic lattice with lattice vectors "e^ = 
(0,0,0),(±1,0,0),(0,±1,0),(0,0,±1),(±1,±1,±1) 
is chosen. Physical variable are defined as moments of 
the distribution functions by 



5> 



pu a 



^2fie ia , Q = ^G 2 . (19) 



The distribution functions evolve in a time step At ac- 
cording to 

At 



— [C fi (x, t, {/<}) + C fi {x + gJAt, t + At, {/*})] , 



(20) 



G i (x + e l At,t + At) - Gi(x,t) 
At 



(16) 



— [C Gi (x,t, {G i })+C Gi (x + e i At,t + At, {G*})] . 

(21) 

This represents free streaming with velocity e*j and a col- 
lision step which allows the distribution to relax towards 
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equilibrium. /* and G* are first order approximations to 
fi(x + tfj At, t + At) and Gj (x + «fj At, t + At) respectively. 
They are obtained by using only the collision operator 
Cfi(x, t, {/»}) on the right of Equation (|20|) and a simi- 
lar substitution in (|2f \ . Discretizing in this way, which 
is similar to a predictor-corrector scheme, has the advan- 
tages that lattice viscosity terms are eliminated to second 
order and that the stability of the scheme is significantly 
improved. 

The collision operators are taken to have the form of a 
single relaxation time Boltzmann eauation|27|. together 
with a forcing term 



cUlt 



Pi 



^ 



P2 



•r 



dx2 



(a) 



— (/i(f,t)-/f(f,t, {/.})) 



C f i(x,t,{fi}) = 
C Gi (x,t,{Gi}) = - — {Gi(x,t)-G?(x,t,{Gi})) 



+Pi(x,t,{fi}), 



(22) 



+M l (x,t,{G l }). 



(23) 



The form of the equations of motion and thermody- 
namic equilibrium follow from the choice of the moments 
of the equilibrium distributions f^ q and G^ q and the 
driving terms pi and Mj. Full details of the algorithm 
can be found in j2jj |^| for two-dimensional and three- 
dimensional simulations respectively. 
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B. Description of the Nanowire 

When a particle is suspended in a nematic liquid, an- 
choring boundary conditions at the particle's surface dic- 
tate the direction of the director field in the vicinity of the 
particle. Using a pinning potential of the form of Eqn iJJJ, 
with n a along the long axis of the wire, ensures the direc- 
tor of the nearby nematic is aligned to the nanowire. The 
nanowire is ferromagnetic, enabling an external magnetic 
field to rotate the wire within the system. 

A numerical algorithm for simulating the inclusion of 
a nanowire in a liquid crystal system is presented in this 
section. We could use an adaptive mesh with boundaries 
aligned on the wire but since the wire diameter is quite 
small (~ 100 nm), we decide against this method due 
to the small resolutions required. For instance, if there 
were 10 mesh points spanning the diameter of wire we 
would be approaching molecular resolution. The alter- 
native is to pixelate a line onto a grid, a problem com- 
mon in computer graphics. First we will start off with 
a two-dimensional description, then expand to three di- 
mensions. We will describe both a standard Brescnham 
line algorithm j3j|, expanded upon by Wu j3||, to pro- 
duce lines on a two-dimensional surface and then we will 
present our current algorithm. 



FIG. 2: (a) Example of the Brensenham anti-aliased line 
drawing algorithm. If the slope is less one, the nodes are posi- 
tioned along vertical lattice intersections, shown as the open 
dots at positions PI and P2. Value of the lattice nodes are 
inversely proportional to the distance from the lattice nodes 
to the wire nodes, (e.g. PI = dx2/dx) For a slope greater 
than one, the horizontal lattice intersections are used, (b) 
Sample cell of a two-dimensional system with an intersecting 
wire. Along the wire, nodes are positioned at fixed distances d 
apart to interact with the lattice nodes. Each node is affected 
by the wire proportionally by the opposite enclosed area rel- 
ative to the cell area. (e.g. P2 is affected from the wire by 
an amount A2/dx 2 .) 



1. Computer Graphics Method 

The standard algorithm for pixelating a line is based 
on Bresenham's algorithm, generalized by Wu [3f| . This 
is illustrated in Fig|2^a). When the slope of the line 
is less than 1, the points where the wire crosses vertical 
grid lines are marked. The value of a particular wire node 
is proportional to the distance of the lattice crossing to 
the opposite node over the cell height. For example, the 
value at P2 = dxi/dx in Figure Ufa). The cndpoints of 
the wire prove to be more cumbersome, involving extrap- 
olation to the next grid line and weighting the points. If 
the slope exceeds one, we swap the roles of x and y and 
exchange the vertical lattice crossings, with horizontal 
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FIG. 3: Three chiral wires are shown together to simulate 
one larger wire. The white spaces are sample locations of the 
nodes along the smaller wires. To avoid commensurability 
with the lattice, the node distance is irrational. Chirality of 
the wire is also included to avoid commensurability issues. 
Each small wire is coiled by 2n about the central axis from 
end to end. 

lattice crossings. This method is easy to implement in 
two-dimensions. Unfortunately it is not easily general- 
ized for three dimensional systems. Other problems arise 
due to the fact that the number of nodes along the wire 
is not conserved as the wire rotates. These issues will 
described in more detail in the next section. 
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FIG. 4: (a) A schematic of the wire immersed in the nematic 
(black dashed) at equilibrium. The magnetic field, the wire, 
and the director field of the nematic are all aligned along the 
same direction, (b) Introducing a magnetic field at angle <f>, 
the wire reorients to an angle 9 balancing the magnetic and 
elastic torques on the wire. This leads to distortions in the 
nematic around the wire. 



2. Node distribution method 

Placing a wire within the simulation will require pixe- 
lation to map it onto the lattice coordinates. Our method 
is illustrated in Fig^b). We divide the wire into a fixed 
number of nodes N. The solid dots are the location of 
the nodes along the wire separated by a distance d, where 
d is incommensurate with the lattice spacings. Typically 
we have used d = 0.24> o dx, where </>o = 0.618 ■ • ■ is the 
golden mean. The wire nodes are then distributed to the 
corners of the cell in which they reside. To distribute the 
wire points to the closest lattice coordinates, the ratio of 
the opposite area from the node over the cell area will be 
calculated (Ex. PI = Al/dx 2 ). It is easy to see that the 
node weight thus defined varies smoothly from 0, when 
the wire node is at the opposite corner of the cell, to 
1, when the wire node goes through the lattice point in 
question. This is done for all nodes along the nanowire. 
If multiple nodes affect a single lattice coordinate, the 
summation of their distributed values is used. Although 
the description above is for a two dimensional system, it 
is easily generalized to three dimensions where ratios of 
opposite volumes to the cell volume is used instead, and 
is the method employed in our simulations. 

In three dimensions an area of fluid will be displaced 
and affected by the wire. To give the wire a lateral di- 
mension greater than dx, we simulated many smaller wire 
strands in unison for the same effect. A diagram of three 
smaller strands forming a single wire is shown in Figure 
[3J In our simulations we use five strands in unison to 
form the nanowire. The white bands represent where the 
nodes lie. Each strand, from end to end, is coiled about 
the central axis by 2tt. This twist on the strands helps 
make the wire nodes even more incommensurate with the 
underlying lattice, something that turns out to be highly 
desirable for producing a smooth response. 

We want the surface energy of the wire to be indepen- 
dent of the way it is discretized. In other words, indepen- 



dent of the number of strands and nodes used (assuming 
a sufficiently high number) and independent of the ori- 
entation with respect to the lattice (i.e. if we rotate the 
numerical mesh but keep all the fields fixed). This re- 
quires that the integral of H w , from Eqn. (|13ll . over the 
surface of the wire be independent of the orientation of 
the wire with respect to the underlying lattice. This will 
be the case if the discretization of § a s dA is a constant 
independent of orientation and wire discretization, where 
a s is proportional to the surface area per node. To cal- 
culate a s , we have 

(b a s dA = 2irrL = V" a s (24) 

where r and L are the radius and length of the wire re- 
spectively, i w is the index for the wire strand, and i n is 
the index for the nodes on wire i w . As a result, i 
the total number of nodes in all the strands. Solving for 
a s gives 



where d is the node spacing along the strands and N s is 
the number of strands used to create the wire. 



IV. RESULTS 
A. Theory 

With no applied field, the long axis of the wire should 
line up parallel with the director, as shown in Figure 01 
(a). Distortions in the orientation of the wire from this 
position, as in Figure 0] (b) will cost elastic energy. The 
director distorts to satisfy the boundary condition at the 
wire surface by introducing an elastic torque T on the 
system. 
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In the equilibrium case where Q is uniaxial, Q a f3 = 
%q(n a np — \5 a p), the elastic free energy density © can 
be rewritten as 

1 



F e i = -[Xi(V-n) 2 + X 2 (n-Vxn) s 
+K 3 (n x V x n) 2 ] . 



(26) 



For analytic work it is useful to work in the approx- 
imation where the elastic constants are all equal. In 
that case, choosing the director field n to have the form 
n x = cos 9 sin ip, n y = sin 9 simp, and n z — costp then 
((26(1 can take the form 

F el = ix[(V^) 2 + sin 2 ^(V0) 2 

+ 2smipn-(V6xVip)]. (27) 

Minimizing l|27|) with respect to the angles leads to two 
equations of equilibrium |37j , 

V 2 V> - sin tJj cos ip( V6») 2 = (28) 
sin ipV 2 9 + 2 cos ipVip • V6» = Q. (29) 

In equilibrium, the nanowire will be positioned in the 
middle of the sample cell parallel to the top and bottom 
surfaces. Using the coordinate system shown in Figure 
^ ip — ^ in equilibrium. If we assume that rotating the 
wire in the plane parallel to the walls produces no out- 
of-plane deformations, then tp will remain constant at 
7r/2. Under this assumption Eqn. (|28|l and (|29|l reduce 
to simply the Laplacian 



V 2 6> = 0. 



If we further assume that the shape of the wire can be 
approximated as a prolate spheroid, the equations can be 
solved analytically. The appropriate coordinate change 
is then 



r± 



r + + r_ 



s = tan 



- (!) • 



(31) 
(32) 
(33) 
(34) 



where a is the interfocal distance. Setting 9 a to a con- 
stant equal to the wire orientation over the surface (i.e. 
m = (cos 9 n , sin 9 a , 0)), the Laplacian l(30() becomes easily 
solvable 133 with solution 



In 



£+1 



in 



(35) 



where the constant £ defines the prolate spheroidal sur- 
face. From equilibrium, we integrate the elastic free en- 
ergy density 1)27(1 over all space, 

U el = y / dV(V9) 2 



= 2-naK9i 



i — i 



In- 



(36) 



which results in an expression for the potential of a pro- 
late spheroidal inclusion in a nematic. 

In the approximation that the three elastic constants 
are equal, K\ — K 2 = K 3 = K, Brochard and de Gennes 
calculated the energy involved in considering a long 
inclusion, such as a wire, within a nematic as a function 
of the inclusions orientation. In the case the anchoring 
at the surface was along the long axis of the wire, they 
exploited an analogy from electrostatics of an object at 
fixed potentials to predict that the energy varies with the 
angle of orientation as 



U ei = 2irCK9 2 , 



(37) 



where C is the capacitance of the wire. Equating l)36|) 
and (|37|l . and making use of Eqn. I|31|) and l|32|) . the 

capacitance is 



C = 



In 



£q + 1 



= 2L\ 1 



In 



(38) 

i 

(39) 



where R and L are the radius and length of the wire 
respectively. In the limit ^ — » 1 the capacitance is simply 
C — > R, that of a sphere, as expected. As j- ~ > 0, the 



(30) capacitance is 



C 



-(f) 



(40) 



which is not the result Brochard and deGennes had at- 
tained, their result containing two incorrect factors of 2. 

Although Eqn. I|40|) is the capacitance for a prolate 
ellipsoid, the wire in our system has a shape of a right 
cylinder. It is not possible to solve Laplace's equation 
exactly for this shape, but expanding in 1/C = l/ln(^), 
Jackson[3^ derived an approximate expression, 



C 



1 1 

c + ? 

1 



l + (l-ln2) 2 - — 
V ' 12 



+ 0(1/C 4 ) .(41) 



For high aspect ratio wires, both Eqn. (|40() and 141(1 
are equivalent. As we shall see below, for low aspect 
ratio wires there is significant differences between the two 
equations and Eqn. 1(41(1 should be used instead. 



B. Torque on a nanowire in the simulations 

In our simulations we will first try to replicate as 
closely as possible the theoretical results of the previ- 
ous subsection. To do this we will first work in the one 
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elastic constant limit and then move on to more realistic 
situations. This will demonstrate the ability of the node 
distribution algorithm to reproduce both the analytic re- 
sults and the results from experiments. 

We start with a wire aligned with the nematic director 
field in equilibrium, as shown in Figure |4*fa). A mag- 
netic field aligned in the direction of the long axis of the 
wire is introduced and slowly rotated in the xy— plane 
(See the schematic in Fig. B£b) and the 3D simulation 
in Fig. n for reference). Taking L 2 = L 3 = 0, to work 
in the one-elastic constant approximation, and if there is 
no out-of-plane distortion, a measurement of the torque 
on the wire from the nematic as a function of rotation 
angle should yield the elastic constant, assuming the ca- 
pacitance of the wire is known. Using Eqn (|37|l for the 
potential energy results in a torque on a wire not aligned 
with the nematic, 



dU 
~d9 



4ttCK9. 



(42) 



Eliminating the dependence on 9, one can simply mea- 
sure the value of the elastic constant via, 



or 



(43) 



Thus, if K is known one can determine C or vice versa. 

Using the two methods described in Sections IIII B II 
and IIII B 21 we can simulate a wire in our system and 
measure the torque imposed by the nematic. The results 
are shown in Fig. [3J It should be immediately appar- 
ent that using the algorithm borrowed from the graphics 
industry, Wu's anti-aliased line method, produces poor 
results (Fig. 0(a)). Aside from being rapidly fluctuat- 
ing, there are profound commensurability effects when 
the wire is at 7r/4 and tt/2. The node distribution al- 
gorithm described in Section ITlI B 21 produces a smooth 
linear relationship between T z and 9, as expected in the 
one-elastic constant approximation |'2l| . as shown in Fig 
EJb). As mentioned in Sections IIII B II Wu's method 
does not conserve the number of nodes along the wire 
(i.e. the wire "density" is not constant). Therefore, as 
the wire rotates the torque, which involves an integral 
over the surface of the wire, fluctuates. In contrast, the 
node distribution algorithm we introduced, conserves the 
wire "density" at all times, as described by Eqn i)24|l . The 
number of nodes along the wire is constant for the algo- 
rithm we proposed, resulting in almost no observable ef- 
fects due to the orientation of the wire with respect to the 
underlying lattice. This leads to a smooth dependence of 
the torque on the angle of the wire with respect to the 
background nematic. Therefore, only the node distribu- 
tion algorithm described in Section IIII B 21 will be used 
for subsequent simulations. 

We are now in a position to test Eqn. |@3J. Within 
the single elastic constant approximation the relationship 
between the capacitance of the wire and the elastic con- 
stant of the nematic is approximately linear, as suggested 
by Eqn 143|) . Using a wire of length 625nm and radius 
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FIG. 5: (a)Wu's algorithm only works in 2-dimensions. The 
system is 1.25/ira x 1.25/im with a plate of length 0.675/im. 
The torque clearly exhibits commensurabilty problems along 
certain orientations of the wire, (b) Using the algorithm 
described in Section IIII B 21 one notices a smooth torque 
along a rotation. The are no indications of any commen- 
surabilty problems as in Wu's algorithm. The system is 
1.25/im x 1.25/im x 1.25/im with a wire of length 0.675/tm,. 
The one elastic constant approximation is used with K = 
10.72piV. To measure the slope a least squares fit is per- 
formed from the to ~2 radians. 



88nm immersed within a nematic in a box of dimensions 
1.25/im x 1.25/im x 1.25/im, the elastic constant of the 
nematic K was varied from 5.2pN to 16.3pA^ in 1.8pN 
increments. For each elastic constant, the wire is rotated 
from to 7r and the slope of the torque versus 9 curve 
is measured (e.g. the slope of the solid line in Fig. |5f b) 
provides one of the points in Fig. Efa)). The resulting 
points are plotted in Figure (a) . As expected, we see 
a linear relationship (the solid line) between K and the 
measured slope of the torque T (diamonds) as expected. 
The slope of this line yields AirC, as can be seen from 
Eqn. g3J. 

Having verified the linear relationship in Eqn. 143|) . 
we can now turn the situation around and use the known 
value of K put into the simulation to measure the un- 
known C . This quantity is not experimentally assessable 
unless an independent measure of the elastic constant K 
is available. We simulated 15 systems where the length 
of the wire was varied from 625nm to 1.875/im and the 
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FIG. 6: (a) The slope of the torque of a 625nm wire as a 
function of the elastic constant. The one-elastic constant ap- 
proximation is used in these simulations, (b) The relationship 
of capacitance as a function of lengths and radii are explored. 
Lengths vary from 625nm to 1.25pim and radii from 44ram to 
198nm. The diamonds are the results from our simulations. 
The exact solution Eq. 1391 and the first order approxima- 
tion, Eq. 1 11)1 for the capacitance of a prolate spheroid are 
the solid and dotted lines respectively. The dashed line is 
the solution of Brochard and deGennes Q] . The dash-dot line 
which traverses through the data is the third order approxi- 
mation, Eqn. 14111 . for the capacitance of a prolate spheroid 
from Jackson. 



diameter from 22nm to 197nm. This resulted in a mini- 
mum aspect ratio of ~ 3 to a maximum of ~ 42 for the 
nanowires in our simulations. To avoid finite size effects, 
we kept the simulation box larger than twice the length 
of the wire (boxes larger than this gave comparable re- 
sults) . The one elastic constant approximation was used 
for all simulations, with K = \Q.72pN. Figure El (b) 
shows our results for C/L as a function of l/ln(2L/i?). 
The dashed line is the incorrect prediction from Brochard 
and de Gennes which clearly does not predict the ca- 
pacitance of a prolate spheroid in a correct manner. In 
Section llV Al we derived the correct expression for the ca- 
pacitance of a prolate spheroid, Eqn. as represented 
by the solid line. The first order approximation is shown 
by the dotted line. Both of these forms for the capac- 
itance converge for high aspect ratio wires, but fail for 
low aspect ratio wires. Interestingly, the first order ap- 



proximation is nearer the simulation results for the wires 
than the exact solution for the capacitance of a prolate 
spheroid. For low aspect ratio wires, the expansion of 
Jackson for a right circular cylinder must be used, and 
one can see close agreement between this third order ap- 
proximation and our results. All three theoretical lines 
(Ecm.l|39 |l -I|41 |l ) converge for high aspect ratio wires. For 
aspect ratios less than ~ 75 Eqn. I|41|) separates from 
both Eqn. I|39|l and Eqn. i|40[l and is obviously the cor- 
rect expression to use. If one were to use the wires as 
micro-rheological probes, they would probably need to 
be a short, low aspect ratio wire (to fit in an in-vivo en- 
vironment) and similar to the sizes studied within this 
paper. 

In real liquid crystals the elastic constants are not equal 
(e.g . for 5CB, K\ and K3 differ by a factor of ~ 1.6 
H3). However, the torque is still observed to be a lin- 
ear function of 8 in experiments [2ll 122^ . This brings up 
the question of which elastic constant, or combination of 
elastic constants, is measured in these experiments. As 
the wire rotates through the nematic, one of the modes 
of deformation, splaying, twisting, or bending, should be 
favorable for the nematic. The distortion field of the ne- 
matic will be greatly affected by the relative differences 
of the strengths in the various elastic constants. Figure 
0is a density plot of the nematic angle with the director 
field of the nematic superimposed for three systems with 
varying elastic constants (A cross section through the 
center of the cell is shown). In all three panels, the wire 
was immersed in the center of the nematic as in figure ^ 
and rotated ^ from equilibrium. In (a) the one elastic 
constant approximation is used with K = 10. 72pN . The 
nematic angle forms homeotropic coronas in the region 
immediately surrounding the wire. In (b) and (c), the one 
elastic constant approximation is no longer valid as var- 
ious elastic constants are used. For (b), K\ = 15.34piV, 
K2 = 3.88piV, and K3 = 8.5pN and in (c) the values 
of K\ and K3 are interchanged. Looking at the density 
plot of the inplane angle, it is discernible there arc clear 
differences in the distortion pattern the nematic uses to 
minimize its free energy for these two systems. The splay 
elastic constant is larger in (b), which explains the dis- 
tortion field of the nematic exhibiting a higher tendency 
to bend within the system. Interchanging the elastic con- 
stants of splay and bend, as in (c), the distortion field is 
dominated by splay, as expected. 

Observations of the in-plane nematic angle show dis- 
cernible differences for the two systems shown Figure 
(b) and (c). However, the torques of the wires rotating 
through the two systems are nearly identical, as shown in 
Figure [S] Even though the elastic constants of splay and 
bend are interchanged, the torques are approximately 
equal implying that for Eqn. I|43|l . the same value of 
K is being measured. Using a least squares fit to the 
data, the solid line in Figure|Sl the slope is ~ 8pN which 
is approximately the value of the smaller of either splay 
or bend depending on the system being observed. Thus, 
experiments will measure the smaller of the splay or bend 
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FIG. 7: Cross sections of three systems of size 3.75/xm x 
3.75[im x 3.75pim with a wire of length 1.875/xm immersed 
in the center. The wire was rotated from equilibrium to an 
angle of ^ . The shading indicates the inplane nematic angle 
with the director field of the nematic superimposed, (a) The 
one elastic constant approximation where Ki = K2 = Kz- 
The director field of the nematic is approximately isotropic in 
the region surrounding the wire. In (b) K\ < Kz and in (c) 
Ki > K 3 




FIG. 8: Comparison of the torques of two systems of size 
3.75/im x 3.75/^m x 3.75/im with a wire of length 1.875/^m. 
The elastic constants of torque represented by the open dia- 
monds are Kx = 15.34piV, K 2 = 3.88piV, and Kz = 8.5pN 
and of the filled triangles are the same elastic constants afore- 
mentioned with K\ and Kz interchanged. The solid line is a 
least squares fit to the data. 



elastic constant and not some average of the elastic con- 
stants as was supposed in Ref. . Of course, it may not 
be known a priori in an experiment which elastic con- 
stant is smaller. However, it should be straightforward 
to observe the birefringence pattern through crossed po- 
larizers and ascertain whether it is more similar to the 
distortion field shown Fig. |7{b) or (c) and thereby ascer- 
tain whether the bend or splay dominates. 



V. CONCLUSIONS 

In this paper we present an algorithm to model a mov- 
ing inclusion in a nematic liquid crystal system. From 
a knowledge of the capacitance of the ferromagnetic 
nanowire, it is possible to deduce elastic constants of the 
liquid crystal from a measurement of the torque. Previ- 
ously quoted expressions for the capacitance of the wire 
Q where found to be incorrect and the correct expres- 
sions, both for a prolate spheroid and for a right circular 
cylinder, were given. Most wires used to date have had 
a very high aspect ratio (~ 100) j^. However, if one 
were to use the wires as an in-vivo probe of the local 
elastic properties in, say, a biological environment then 
much lower aspect wires, or prolate spheroids may be 
more appropriate. In that case, (i.e. for aspect ratios 
less than around 75) either the exact expression for the 
prolate spheroid or the third order expansion by Jackson 
for the appropriate capacitance should be used. 

For realistic systems with multiple elastic constants, we 
found that the system will tend to trade between splay 
or bend distortions depending on which one cost less free 
energy. As a result, the linear elastic response of the 
nematic to a rotation of the wire is dominated by the 
smaller of the splay or bend elastic constants. Thus, 
experiments will typically measure the smaller of K\ or 
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K$ and not some combination. It should be possible to 
deduce the elastic constant being measured based on the 
distortion field of the nematic around the rotated wire, 
as demonstrated in Fig. [7| 

In future works, we will examine the nonlinear response 
that occurs at wire rotations close to n and the dynamical 
response. 
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